% This code can be executed to generate Figure 7. Similar to the code for
% Figure 5A, The purple curve is theoretical peeling maximum (PeelF) is 
% expressed as a function of the junction angle. The yellow curve is the 
% elastic force (intersectUFelas) where the peeling and elastic energies 
% are equal (intersectU). The orange curve is point where the elastic model 
% becomes linear at x=alpha (linearpts). The tensile strength of the tape 
% (TS_tape) is the black dashed line. The grey histogram (hist) shows the
% measured junction angles of the natrual recluse spider's loops. The
% vertical black line (avgangle) is the average of these values. The
% vertical blue dashed lines mark the boundaries between the failure
% regions. Additional formatting is done externally to exactly match the 
% published Figure 7.
% The folder titled, "Raw Data Files," must be open to run this code!
clear;

%Constants for Recluse Silk
t=50E-3; %um
w=7; %um
E=15E-3; %N/um2
gamma=47/1E9;%N/um
beta=200; %200 loops/cm = 50 um/loop
TS_tape=0.629E-3; 
std_TS=0.1895E-3;
slope=(E*w)./1000; %N/mm

%Determine Energy Equivalence and Transition to Linearity (alpha)
angles=linspace(1.25,76,10000);
intersectUFelas=zeros(1,length(angles));
xmaxs=zeros(1,length(angles));
TS_tapeFelasx=zeros(1,length(angles));
linearpts=zeros(1,length(angles));
linearptsFelas=zeros(1,length(angles));

for i=1:length(angles)
    phi=angles(i);
    xmax=2.*w.*tand(phi./2);
    xmaxs(i)=xmax;
    %Peeling Energy
    Upeelinc=@(x) 2.*gamma.*((cscd(phi./2)).^2).*(1./2).*(x.^2).*(cotd(phi./2));
    Upeeldec=@(x) 2.*gamma.*((cscd(phi./2)).^2).*((1)./2).*((-1.*(xmax-x).^2)+2.*(xmax./2).^2).*(cotd(phi./2));
    %determine peeling energy for increaing and decreasing sections of
    %force curve
    xx=linspace(0,20,10000);
    for j=1:length(xx)
        pos=xx(j);
        if pos>=(xmax/2)
            xxinc=xx(1:j);
            Upeelincvec=Upeelinc(xxinc);
            half=j;
            break
        end
    end
    for h=half:length(xx)
        pos=xx(h);
        if pos>=xmax
            xxdec=xx(half:h);
            Upeeldecvec=Upeeldec(xxdec);
            break
        end
    end
    %add tail of constant value to make sure vectors are same length
    tailU=ones(1,(length(xx)-h)-1).*max(Upeeldecvec);
    Upeel=[Upeelincvec Upeeldecvec tailU];
    
    %Elastic Force
    Felas=@(x) -(1/2).*E.*t.*x.*tand(phi).*log(abs((2.*x.*cotd(phi./2)-beta.*tand(phi))./(-beta.*tand(phi))));
    Felas_val=Felas(xx);
    for n=1:(length(Felas_val)-1)
        m=(Felas_val(n+1)-Felas_val(n))/(xx(n+1)-xx(n));
        b=Felas_val(n)-m*xx(n);
        %linearize elastic force
        if m>slope
            xremain=linspace(xx(n),max(xx),(length(xx)-n+1));
            Felas_val(n:length(xx))=xremain.*slope+b;
            break
        end
        linearpts(i)=xx(n);
    end
    %Elastic Energy
    Uelas=cumtrapz(xx,Felas_val);
    
    %Difference between peeling and elastic energies to determine where the
    %peeling and elastic energies intersect
    subU=Uelas-Upeel;
    for k=1:length(subU)-1
        test=subU(k);
        if test>0
            intersectUFelas(i)=Felas_val(k);
            break
        end
    end

    %Determine Alpha
    linearptsFelas(i)=Felas_val(n);
end

%Determine Peeling Model
PeelF=@(phi) 2.*gamma.*w.*((cscd(phi./2)).^2);
phiphi=linspace(0.6,89,1000);
peel_f=PeelF(angles)./areacalculation(w,angles);

%Determine Boundary Between Regions II and III and Regions I and II
linearptsFelas=linearptsFelas./areacalculation(w,angles);
intersectUFelas=intersectUFelas./areacalculation(w,angles);
for i=1:length(angles)
    alphaloc=linearptsFelas(i);
    peelloc=peel_f(i);
    if peelloc<alphaloc
        R_I_II=(angles(i)+angles(i-1))/2;
        break
    end
end
for i=1:length(angles)
    peelloc=peel_f(i);
    eqUloc=intersectUFelas(i);
    if angles(i)>R_I_II && peelloc>eqUloc
        R_II_III=(angles(i)+angles(i-1))/2;
        break
    end
end
for i=1:length(angles)
    alphaloc=linearptsFelas(i);
    eqUloc=intersectUFelas(i);
    if angles(i)>R_I_II && eqUloc<alphaloc
        R_II=(angles(i)+angles(i-1))/2;
        break 
    end
end

%Determine the meansured loop junction angles of natural recluse spider
%loops
Data=importdata('Loxosceles Loop Junction Angle Measurements.txt');
angleslox=Data.data;
avgangle=mean(angleslox);
SD=std(angleslox);

figure;
hold on;
yyaxis right
ylim([0 0.08]);

%Plot average junction angle recluse spiders make their loops
avglin=ones(1,1000).*avgangle;
avglinvec=linspace(0,50,1000);
plot(avglin,avglinvec,'k','Marker','none','LineWidth',2);

%Plot histogram of all loop junction angles measured
%Note: This is the same histogram in Figure 2C
hist1=histogram(angleslox,'BinWidth',5,'Normalization','pdf',...
    'FaceColor',[150/255 150/255 150/255],'FaceAlpha',1/3,'EdgeColor',[0 0 0],'EdgeAlpha',1);

%Optional: plot plus/minus one standard deviation of the mean junction
%angle that recluse spiders make their loops
%patch([avgangle-SD avgangle-SD avgangle+SD avgangle+SD],[0 20 20 0],...
%    'g','FaceAlpha',0.2);

ylabel('Counts','FontSize',14);
yyaxis left

%Plot Peeling Model
plot(phiphi,PeelF(phiphi)./areacalculation(w,phiphi),'Marker','none','Color',[127/255 0 255/255],'LineStyle','-','LineWidth',2);

%Plot Tensile Strength of Tape
hline=ones(1,length(phiphi)).*TS_tape;
plot(phiphi,hline./areacalculation(w,phiphi),'Color','k','LineStyle','--');

%Plot where Elastic Energy=Peeling Energy
plot(angles,intersectUFelas,'Color',[255/255 193/255 37/255],'LineWidth',2,'LineStyle','-');

%Plot where slope of elastic force=modulus of single tape (alpha)
plot(angles,linearptsFelas,'Color',[255/255 126/255 24/255],'LineWidth',2,'LineStyle','-');

%Determine Extrapolated Range of allowed junction strengths
n_iterations=1;
norm_LOE_extrap=zeros(1,n_iterations);
for l=1:n_iterations
    %Determine average Loop Strength based on large # distribution of Loxo Hist
    %determine probability density function (pdf) of measured loxo angles
    f=@(sigma,mu,alpha,x) ((2./(sigma.*sqrt(2*pi))).*exp(-((x-mu).^2)./(2.*sigma.^2))).*((1/2).*(1+erf((alpha.*((x-mu)./sigma))./sqrt(2))));
    sigma=SD;
    mu=avgangle;
    alpha=2.5;
    test=hist1.BinEdges;
    valuesx=zeros(1,(length(test)-1));
    for k=1:length(valuesx)
        valuesx(k)=(test(k)+test(k+1))/2;
    end
    valuesy=hist1.Values;
    ecounts=ones(length(valuesy),1).*1;
    [fitresult]=fit(valuesx',valuesy',f,'StartPoint',[sigma,mu,alpha],'Weights',power(ecounts,-2.0));
    pdf_angleslox=f(fitresult.sigma,fitresult.mu,fitresult.alpha,angles);
    %Optional: plot pdf of measured loxo angles
    %yyaxis right
    %plot(angles,pdf_angleslox,'-r');
    %yyaxis left
    
    %Determine number of loops at a particular junction angle based on
    %multiplying a N_loops with the pdf
    theo_lox_angle_dist=pdf_angleslox.*8.108; %N_loops=10,000    8.108=N_loops=1,000
    %assemble a vector of the theoretical loop junction angles
    %(theo_lox_angles)
    bins=linspace(min(angles),max(angles),101);
    number_of_junctions=zeros(1,length(bins)-1);
    theo_lox_angles=[];
    for i=1:(length(bins)-1)
        left_edge=bins(i);
        right_edge=bins(i+1);
        for j=1:length(theo_lox_angle_dist)
            if angles(j)>=left_edge && angles(j)<right_edge
                number_of_junctions(i)=number_of_junctions(i)+theo_lox_angle_dist(j);
            end
        end
        for k=1:round(number_of_junctions(i))
            theo_lox_angles=[theo_lox_angles left_edge+(right_edge-left_edge)*rand();];
        end
    end
    
    %Determine theoretical junction strength of theoretical distribution of
    %loop angles
    theo_lox_TS_l=[];
    theo_lox_TS_n=[];
    theo_lox_TS=[];
    theo_lox_TS_langles=[];
    theo_lox_TS_nangles=[];
    theo_lox_TS_angles=[];
    n_l_prob=0.15; %probability to be in nestled conformation
    for i=1:length(theo_lox_angles)
        angle=theo_lox_angles(i);
        n_or_l=rand(); 
        if angle<=R_I_II
            upper=(TS_tape+std_TS)./areacalculation(w,angle);
            lower=(TS_tape-std_TS)./areacalculation(w,angle);
            theo_lox_TS=[theo_lox_TS lower+(upper-lower)*rand()];
            theo_lox_TS_angles=[theo_lox_TS_angles angle];
        elseif angle>R_I_II && angle<=R_II_III && n_or_l>n_l_prob
            upper=(TS_tape+std_TS)./areacalculation(w,angle);
            for m=1:length(angles)
                if angle<=angles(m)
                    break
                end
            end
            lower=linearptsFelas(m);
            theo_lox_TS_l=[theo_lox_TS_l lower+(upper-lower)*rand()];
            theo_lox_TS_langles=[theo_lox_TS_langles angle];
        elseif angle>R_II_III && n_or_l>n_l_prob
            upper=7*PeelF(angle)./areacalculation(w,angle);
            lower=PeelF(angle)./areacalculation(w,angle);
            theo_lox_TS_l=[theo_lox_TS_l lower+(upper-lower)*rand()];
            theo_lox_TS_langles=[theo_lox_TS_langles angle];
        elseif angle>R_I_II && angle<=R_II
            for m=1:length(angles)
                if angle<=angles(m)
                    break
                end
            end
            upper=linearptsFelas(m);
            lower=PeelF(angle)./areacalculation(w,angle)-0.25*PeelF(angle)./areacalculation(w,angle);
            theo_lox_TS_n=[theo_lox_TS_n lower+(upper-lower)*rand()];
            theo_lox_TS_nangles=[theo_lox_TS_nangles angle];
        elseif angle>R_II && angle<=R_II_III
            for m=1:length(angles)
                if angle<=angles(m)
                    break
                end
            end
            upper=intersectUFelas(m);
            lower=PeelF(angle)./areacalculation(w,angle)-0.25*PeelF(angle)./areacalculation(w,angle);
            theo_lox_TS_n=[theo_lox_TS_n lower+(upper-lower)*rand()];
            theo_lox_TS_nangles=[theo_lox_TS_nangles angle];
        elseif angle>R_II_III
            upper=PeelF(angle)./areacalculation(w,angle)+0.25*PeelF(angle)./areacalculation(w,angle);
            lower=PeelF(angle)./areacalculation(w,angle)-0.25*PeelF(angle)./areacalculation(w,angle);
            theo_lox_TS_n=[theo_lox_TS_n lower+(upper-lower)*rand()];
            theo_lox_TS_nangles=[theo_lox_TS_nangles angle];
        end
    end
    %scatter these modeled data points and color according to Table 1
    scatter(theo_lox_TS_angles,theo_lox_TS,30,'filled','MarkerFaceAlpha',1/3,'MarkerFaceColor',[0 0 0]...
            ,'MarkerEdgeColor',[0 0 0],'MarkerEdgeAlpha',2/3);
    scatter(theo_lox_TS_langles,theo_lox_TS_l,30,'filled','MarkerFaceAlpha',1/3,'MarkerFaceColor',[255/255 112/255 101/255]...
            ,'MarkerEdgeColor',[0 0 0],'MarkerEdgeAlpha',2/3);
    scatter(theo_lox_TS_nangles,theo_lox_TS_n,30,'filled','MarkerFaceAlpha',1/3,'MarkerFaceColor',[0/255 190/255 30/255]...
            ,'MarkerEdgeColor',[0/255 100/255 0/255],'MarkerEdgeAlpha',2/3);
    %Determine the modeled junction strength
    all_theo_lox_TS=mean([theo_lox_TS theo_lox_TS_l theo_lox_TS_n]);
    norm_LOE_extrap(l)=mean([theo_lox_TS./(TS_tape./areacalculation(w,theo_lox_TS_angles))...
        theo_lox_TS_l./(TS_tape./areacalculation(w,theo_lox_TS_langles))...
        theo_lox_TS_n./(TS_tape./areacalculation(w,theo_lox_TS_nangles))]);
end
%
%Plot Region Boundaries
line1_2=ones(1,1000).*R_I_II;
line2_3=ones(1,1000).*R_II_III;
yy=linspace(0,0.1,1000);
plot(line1_2,yy,'b','LineStyle','--','Marker','none');
plot(line2_3,yy,'b','LineStyle','--','Marker','none');

%Formatting the plot
xlabel('Junction Angle (\circ)','FontSize',14); 
ylabel('Max Load/Junction Area (MPa)','FontSize',14);
ax = gca;
ax.XAxis.MinorTick = 'on';
ax.XAxis.MinorTickValues = ax.XAxis.Limits(1):2.5:ax.XAxis.Limits(2);
ax.FontSize=14;
ylim([0 0.7E-5]);
xlim([0 90]);
hold off;
%%
%Plot distribution of junction strengths similar to Figure 2B
R1_norm_LS=theo_lox_TS./(TS_tape./areacalculation(w,theo_lox_TS_angles));
l_norm_LS=theo_lox_TS_l./(TS_tape./areacalculation(w,theo_lox_TS_langles));
n_norm_LS=theo_lox_TS_n./(TS_tape./areacalculation(w,theo_lox_TS_nangles));
norm_LS=[R1_norm_LS l_norm_LS n_norm_LS];
figure;
hold on
xx=ones(1,length(norm_LS)).*(mean(norm_LS).*100);
yy=linspace(0,600,length(norm_LS));
plot(xx,yy)
hist3=histogram(norm_LS.*100);